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^ • Abstract: 

■ Motivated by the sign problem in several systems, we have developed a geometric simulation 
V/-^ . algorithm, based on the strong coupling expansion, which can be applied to abelian pure 

, gauge models. We have studied the algorithm in the U{1) model in 3 and 4 dimensions, 

I and seen that it is practical and is similarly efficient to the standard heat-bath algorithm, 
but without the ergodicity problems which comes from the presence of vortices. We have 

^ ^ also applied the algorithm to the Ising gauge model at the critical point, and we find hints 

■ of a better asymptotic behaviour of the autocorrelation time, which therefore suggests the 



possibility of a smaller dynamical critical exponent with respect to the standard heat-bath 
algorithm. 
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1. Introduction 

The expected scenario for the phase diagram of QCD in the chemical potential-temperatu- 
re plane is based on the conjecture that, besides the well established hadronic and quark- 
gluon plasma phases, there exists a new state of matter. This new phase is characteristic 
of the high density, low temperature regime: the asymptotic freedom of QCD and the 
known instability of large Fermi spheres, in the presence of weak attractive forces, results 
in a pairing of quarks with momenta near the Fermi surface (in analogy with the Cooper 
pairing in solid state systems at low temperature). A condensation of quark pairs should 
be the distinctive signal of the new phase, and it has indeed been predicted using simplified 
phenomenological models of QCD. Unfortunately the lattice approach, the most powerful 
tool to perform first principles, non perturbative studies, is afflicted in the case of finite 
density QCD by the well known sign problem that, notwithstanding a large amount of work 
by many people [1], has prevented until now any significant step towards the understanding 
of this new phase. 

An alternative strategy to overcome the sign problem consists in reformulating the 
theory, not in terms of the microscopic degrees of freedom (quarks and gluons in the QCD 
case), but in perhaps some other more physical variables (mesons, baryons, etc. [2, 3]). 
Following this line, Karsch and Miitter [2] analyzed the phase diagram of strongly coupled 
QCD with SU(3) gauge fields, staggered fermions and four flavors. By integrating out the 
SU(3) gauge fields [4], Karsch and Miitter constructed a new representation of the QCD 
partition function as a statistical system of monomers, dimers and polymers (MDP-system). 
Even if at finite density some of the Boltzmann weights were still negative in the new QCD 
representation, they claimed that the dominant contributions to the partition function 
had positive weights; and this allowed for the construction of an algorithm that generated 
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configurations distributed according to the absolute value of the integration measure, with 
the sign of the Boltzmann weights being absorbed into the observables [5, 6]. 

Unfortunately the continuum limit of QCD lies not in the strong coupling regime but 
in the weak coupling regime, and here the pure gauge Boltzmann factor strongly changes 
the gauge integration rules [4] , invalidating the MDP-QCD representation for the partition 
function. Indeed the MDP representation would be the zeroth order expansion of the 
partition function in powers of /3, the square of the inverse gauge coupling. 

The Karsch and Miitter results for the strong coupling limit of finite density QCD 
could, in principle, be extended to higher order in the expansion in (3. This extension will 
obviously increase the complexity of the graphs contributing to the new representation 
of the QCD partition function, and fTirtlicrmore, it is not clear at all whether the sign 
problem in the new representation will become a weak sign problem, as in the strong 
coupling case, or a severe sign problem [7]. In any case, the physical relevance of the 
system under consideration makes it worthwhile to try such a line of research. In addition 
there are other physically relevant gauge models with a sign problem, such as systems with 
a topological O-teim in the action, where these techniques could also be applied [8-17]. 

Several years ago Prokofev and Svistunov [18] proposed new algorithms for classical 
statistical spin systems, the worm algorithms, based on the use of a new representation of 
the partition function for these systems, which emerges from the high temperature expan- 
sion. The use of this expansion involves introducing a new configuration space for the Ising 
model made up from monomers, dimers and closed paths, with their corresponding exclu- 
sion rules. Prokofev and Svistunov showed how these worm algorithms applied to several 
spin systems essentially eliminate critical slowing down, with dynamic critical exponents 
close to zero, and yet remain local. These results have been recently corroborated in [19] 
for the Ising model , reporting a detailed error analysis and a generalization of the method 
of [18]. 

Our aim in this paper is to construct and apply algorithms for the strong coupling 
representation of the partition function of three and four-dimensional abelian pure gauge 
models, in particular U{1) and Z{2) [20]. Independently of our aforementioned motivation, 
it seems interesting to us to construct these algorithms and to compare them with standard 
ones, such as heat-bath, in order to see if, as reported in [18, 19, 21] for spin systems, critical 
slowing down is also practically absent in gauge systems. This is precisely the reason why we 
have also simulated the three-dimensional Ising-gauge model. Indeed, contrary to the U{1) 
model, the three-dimensional Ising-gauge model, which is dual of the Ising model in three 
dimensions, shows a second order phase transition, were critical slowing down is present 
in the standard algorithms. The generalization of the strong coupling representation of 
the partition function to the U{1) model with dynamical fermion matter fields, which in 
principle suffers from the sign problem (at least in more than two dimensions) [22, 23], will 
be reported in a forthcoming publication. 

The paper is organized as follows. In section II we analyze the strong coupling repre- 
sentation of the partition function, discuss the possible structures which can appear in the 
abelian model, and compute the configuration transition probabilities for local changes. 
Section III shows how to compute thermodynamic quantities in the new representation, 
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such as the plaquette energy, specific heat, Wilson loops and correlation functions. In 
section IV we describe the implementation of our algorithm, report numerical results for 
the three and four dimensional U{1) compact model, and compare the efficiency of our 
algorithm against heat-bath. Section V contains our results for the three-dimensional 
Ising-gauge model which, as stated before, is a good laboratory to check the efficiency 
of the algorithm near second order phase transition points. Last section is devoted to 
summarize our conclusions and discuss possible extensions of this work. 



2. Abelian groups 



We will consider an (euclidean) abelian U{1) gauge theory on the lattice with the usual 
Wilson action, 

S = -pRe J2 UMU.{n + f,)Ulin + i.)C/t(n) = - ^ + [/*) (2.1) 



IJi<V 



where Np is the number of plaquettes of the lattice, k indexes the plaquettes of the lattice 
and Uk and are the oriented product of gauge fields along plaquette k and its complex 
conjugate respectively. 

The starting point of the algorithm is the strong-coupling expansion of the partition 
function: 

z=j [\{duM] n = / [m n (f; ^(c/. + u^t] = 

n,n k=l k ln=0 ' J 
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If we expand the product over the plaquettes, we get a sum of terms, which can be uniquely 
labeled by giving a couple of natural numbers for each plaquette (n^, n^), where and 
correspond to the powers of Up and U*, and the super- index a refers to a specific term in the 
expansion. After integration over the gauge fields, the only non-vanishing contributions to 
Z are those in which every link corresponding to a plaquette with (n^, n^) 7^ (0, 0) appears 
accompanied by its complex conjugate. Each of these contributions therefore can consist 
only of plaquettes with = n' 
can write 



^ , and plaquettes which form closed orientable surfaces. We 
Z{P)=Y,C%P) (2.3) 



where a indexes the set of non-vanishing contributions. The (non-negative) weight corre- 
sponding to such a contribution is given by 



c"(/3) = n 



(2.4) 



k- 
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Now we can view Z as the partition function of a new system, where a configuration a is 
characterized by the set of numbers (n^, n^) corresponding to a non-vanishing contribution 
to the strong couphng expansion, with probability w'^ = C^/Z, and implement a Monte 
Carlo in this new configuration space. We define elementary Monte Carlo moves consisting 
in adding a "double plaquette" , (n^, n^) — (n^ + 1, + 1) (as well as the opposite move, if 
both ra^,^^ > 0), and adding an elementary cube (as well as removing it if present). This 
ensures that we stay inside the set of non- vanishing contributions. We choose the transition 
probabilities to be proportional to the weights of the corresponding configurations in the 
partition function, which ensures detailed balance. This can be done efficiently, and it is 
also easy to see that the algorithm is ergodic. ^ 

The gauge Zp model can be handled in the same way, the only difference being that 
now there are additional non-vanishing contributions given by powers of p of the plaquette 
or its complex conjugate, that is, we must allow contributions with n'^ — n'^ = (mod p). 
The weight of a configuration is given by the same expression as before (2.4). The only 
modification required in the algorithm is adding a new move which creates or destroys a p 
power of the plaquette or its complex conjugate. 



3. Observables 

The computation of observables is quite easy in this representation. The definition of the 
observable plaquette^ is 

{Ih) = ^dpln{Z) (3.1) 



Np 



Using (2.3) and (2.4), we obtain 



\aeC / ^ aeC ^ 



with 



n 



Np 

fc=i 



n 



Np 
k=l 

The quantity C^/Z is the probability u;" of each configuration a. The mean of the ob- 
servable plaquette is then 

(^□) = ;a^E^"(^" + ^") (3-3) 



single surface that wraps around the lattice cannot be created by this algorithm starting from the 

empty configuration, but this is only a finite- volume correction. 

^We must remark that the observable plaquette and the plaquettes living in our lattice are not the same, 
although they are strongly related. The observable plaquette refers to the minimal Wilson loop, whereas 
the plaquettes refers to geometric entities, living on the lattice. In order to keep the discussion clear, we 
will always refer to the minimal Wilson loop as observable plaquette. 
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that is to say, the observable plaquette is equal to the mean value of the sum of the 
occupation numbers n" plus n", divided by /3 and normalized by Np. 
Another interesting observable is the specific heat: 

Cv = dp{Pb) (3.4) 
We can profit from the previous expression of Pn (3.3) to find the following equation 



= i E - (e »"^) ' - E 



(3.5) 



Sometimes, it is interesting to compute correlation observables, such as the Wilson 
loop (larger than the single observable plaquette), or the plaquette-plaquette correlation 
function. To compute these it will prove helpful to introduce a pair of variable coupling 
constants ^k}, which depend on the plaquette site k, in such a way that the partition 
function reads now 

m,Pj) = I mUei^'''^^''*^) (3.6) 

k 

The weight of the configurations changes accordingly 

Np r.nl rfi 

CiPiJj) = TT „^-a^ \ , r (3.7) 

Now the correlation functions or the Wilson loops are computed by simple derivation, and 
then taking all the Pk, Pk to the same value. For instance, the 2x1 Wilson loop can be 
calculated as 

where k and k+1 are contiguous plaquette sites. The generalization of this result to larger 
planar Wilson loops is straightforward. Indeed, the expectation value of any planar Wilson 
loop can be computed as the mean value of the product, for all plaquettcs enclosed by 
the loop, of the occupation number of plaquettes in each plaquette site, multiplied by a 
factor 2//? to a power which is just the number of plaquettes enclosed by the loop. This 
observable is quite remarkable, for it is computed as a product of occupation numbers, 
and features a (2//?)"^''^" factor, which eventually may become exponentially large or small 
as the size of the loop increases. All these particular facts render this observable hard to 
compute, as we will see in the numerical results. 

In the same way we can also obtain the correlation functions for two arbitrary plaque- 
ttes on the lattice: 
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and 

Finally from (3.9) and (3.10), and taking into account the symmetry of the model, we can 
write: 

(ReUkReUi)^ = {{uk + fik) {ni + fii)) - {uk + n^) {ni + ni) (3.11) 



(ImC/fclmC/;)^ = --^ ((nfc - n^) (n; - ni)) (3.12) 



where the brackets denote average over configurations. 



4. Implementation and numerical results 

In order to see our algorithm (which we should call, from now on, geometric algorithm) at 
work, we have performed numerical simulations of several lattice gauge theory systems in 
three and four dimensions. Our aim is to check the goodness of our approach, comparing 
the results we obtain using the geometric algorithm with those coming from more standard 
ones; hence we want to compare the properties of the algorithm itself, in terms of numerical 
efficiency and autocorrelation times with, for example, the usual heat-bath algorithm. We 
have in mind the results of [19], where it was claimed that, with a similar algorithm, at a 
second order critical point in a non-gauge system, there is a very strong reduction in the 
critical slowing down. 

Let us start with the three dimensional U(l) lattice gauge model: this model is known 
to have a single phase from strong to weak coupling. We have chosen to measure two simple 
observables, namely the plaquette observable and the specific heat, following the definitions 
given in the preceeding Section. We have simulated the model with our algorithm and with 
a standard heat-bath for a large interval of (3 values using a 12^ lattice; we allowed the 
system to reach thermalization for 5 x 10^ iterations and then measured the observables for 
10® iterations. Errors were evaluated using a Jackknife procedure. The results are shown 
in Fig 1. We can easily see in this figure that the two simulations give essentially the same 
results. 

Almost the same situation can be depicted also for the four dimensional U(l) model; 
the results of a similar set of simulations, performed with the two algorithms on a 16^ lattice, 
are shown in Fig 2. Here the only difference can be seen near the phase transition point. 
Remember that due to the difference in finite volume terms between the two algorithms, 
the precise pseudo-critical coupling value at finite volume has to be slightly different. We 
have calculated a few more points very close to the critical beta for each algorithm, and 
the results can be seen in tables 1 and 2. Note the much larger value of the specific heat 
for some of these points as compared with fig.2. 

These results allow us to infer that the geometric approach is able to reproduce all the 
features of the models under investigation, and when differences are seen, they can be easily 
explained on the difference between finite volume terms. To study more carefully these 
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Figure 1: Three dimensional U(l) lattice gauge system. Errors are smaller than symbols. 
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Figure 2: Four dimensional U(l) lattice gauge system; please note the different scale for the specific 
heat on the right. Largest errors (those on the pseudo-critical point) are smaller than symbols. 



differences we have calculated PdL), the critical coupling for each algorithm at different 
lattice sizes. We present in table 3 the results. We also include in the table the results of a 
fit of Pc{L) for the three largest lattices from each set to the finite-size scaling law expected 
for a first-order phase transition [24]. This gives a good fit and consistent results for the 
infinite volume limit. 
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Plaquette 



1.010700 
1.010850 
1.010900 
1.011100 
1.011600 



0.62545(34) 
0.64014(143) 
0.64940(54) 

0.65473(6) 
0.65576(12) 



24.1(3) 
88.2(2) 
47.9(4) 

2.6034(1) 
2.5559(2) 



Table 1: Heat-bath results near /3c 



Plaquette 



1.011120 0.6373(45) 67.2(7) 

1.011160 0.6423(46) 67.4(3.6) 
1.011420 0.6549(1) 2.9660(2) 

Table 2: Geometric results near /3c. 
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Figure 3: Real part of the Wilson loop versus the loop area for the confining phase (/3 = 0.9) 



in the four-dimensional U(l) gauge model, 
algorithm. The lattice volume was 16^. 



Notice the absence of the 5x5 loop in the geometric 



The presence of two clearly differ- 
ent phases in this model, namely a con- 
fining and a Coulomb one, allows us to 
study the behaviour of the Wilson loop 
results in two different physical situa- 
tions; as above, we have also performed 
standard simulations for a cross check 
between the two approaches. In Figs. 
3 and 4 we report the behaviour of the 
Wilson loop in both phases (confining 
in Fig. 3 and Coulomb in Fig. 4) and 
in lattices of different size (12^, 14"^, 
16^). 



L Pc{L) (heat-bath) (5c{L) (geometric) 
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1.00171(8) 


1.00878(20) 


10 


1.00936(11) 


1.01062(7) 


12 


1.01027(8) 


1.01100(7) 


14 


1.01064(4) 


1.01103(20) 


16 


1.01084(17) 


1.01116(14) 


oo 


1.01108(10) 


1.01120(22) 



Table 3: (3c{L) for heat-bath and geometric algorithm 
respectively. 
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Figure 4: Real part of Wilson loop versus the loop perimeter for the Coulomb phase (/3 = 1.1) in 
the four-dimensional U(l) gauge model. The lattice volume was 16^. 

These figures deserve some comments. The geometric algorithm seems to suffer from 
larger statistical errors than the heat-bath method, regardless of the phase of the system. 
To understand this result, we should have a close look at the inner machinery of the al- 
gorithm, in particular, the way the Wilson loop is computed (see eq. (3.8)). First of all, 
the mean value of the Wilson loop is computed as a sum of integer products, implying the 
existence of large fluctuations between configurations. For example, doubling the occupa- 
tion number of a single plaquette doubles the value of the loop. This is a quite common 
fiuctuation at the /3 values of our simulations, and the fiuctuations will increase as the loop 
(and therefore the number of plaquettes) grows. To complicate the computation further, 
we are trying to calculate an exponentially small quantity by summing integer numbers. 
The discrete nature of this computation tells us that non-zero values of the quantity must 
appear with an exponentially small probability. This explains the inherent difficulties of 
the large Wilson loops (4x4 and greater) measurement in the confining phase. The result 
is shown in Fig. 3: the mean value of the 5x5 Wilson loop was exactly zero in the geomet- 
ric algorithm, which is of course wrong. Finally, the expectation value of the Wilson loop 
is proportional to a {2/(5)^ factor, with A the loop area. This value may become huge (or 
tiny) for large loops and low (or high) values of beta, thus enhancing the problems that 
arise from the discreteness of the algorithm ^. 

Notwithstanding the stronger fluctuations in the large Wilson loops within the geo- 
metric algorithm discussed above, it has a clear advantage against heat-bath: it does not 
suffer from ergodicity problems. Indeed the results for the Wilson loop at /3 = 3 reported 

^See [25, 26] for a discussion of some numerically efficient algorithms for the calculation of large Wilson 
loops and Polyakov loop correlators. 
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Figure 5: Real part of Wilson loop versus the loop perimeter for a large /3 value {(3 = 3.0) in the 

four-dimcnsional U(l) gauge model. Notice the difference in performance between the hot and the 
cold starts of the heat-bath algorithm. The lattice volume was 12^. 



in Fig. 5 strongly support the previous statement. The points obtained with the geometric 
algorithm nicely follow the weak coupling prediction of [27], whereas the heat-bath results 
for large Wilson loops, obtained from a hot start, clearly deviate from the analytical weak 
coupling prediction. The origin of this anomalous behavior in the heat-bath case is related 
to the formation of vortices, which are metastable states, that become extremely long lived 
in the Coulomb phase [28]. 

We have also calculated the plaquette-plaquette correlation (of the real (3.11) and of 
the imaginary (3.12) parts) in both phases and for plaquettes living in the same plane. 
Here we expect a much milder behaviour for the geometrical algorithm, for there is no 
large (2//3)^ factor, and the fluctuations are reduced to a couple of plaquettes. The results 
are shown in Figs. 6, 7, 8 and 9. 

In all the cases, the numerical results obtained with the geometric and heat-bath algo- 
rithms essentially agree, except for the correlations of the imaginary part of the plaquettes 
in the Coulomb phase (Fig. 9) , where a clear discrepancy for distances larger or equal than 
4 is observed. Again in this case the reason for this discrepancy is related, as for the Wilson 
loop results previously discussed, to the formation of extremely long-lived metastable states 
[28] in the heat-bath simulations, which seem to be absent in the geometric algorithm. In- 
deed we have verified, with simulations in 12^ lattices, that when we start the heat-bath 
runs from a cold configuration, the disagreement on the correlations of the imaginary part 
of the plaquettes in Coulomb phase at large distances basically disappear. There are still 
small discrepancies in this case, but they can be reasonably attributed to the difference in 
finite volume terms between the two algorithms. 
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Figure 6: Correlation function of the real part of the plaquette versus plaquette-plaquette distance 
in lattice units, for the four-dimensional U(l) lattice gauge model in the confining phase (/3 = 0.9). 
Beyond distance 4, the error became far larger than the expectation value of the correlation. The 
lattice volume was 16*. 



To compare computational costs we define a figure of merit, which is the product of the 
squared error times the cpu time. We expect the error to vanish like 1 / V-^Monte Carlo i ^iid 
therefore the quantity defined above should tend asymptotically to a constant. We show 
the value of this quantity for several observables in both phases and for both algorithms 
in Fig. 10. We can see that the performance of both algorithms is quite comparable. 
The differences that are seen could conceivably change if one were to optimize the specific 
implementations, but none is obviously much more efficient than the other for the models 
studied. 

In particular, for the plaquette observable and the specific heat, both algorithms have 
a similar figure of merit. Prom our point of view, the differences are not quite significant, 
and could change with careful optimizations. The real plaquette-plaquette correlation is 
quite another story, for the differences become significative in the Coulomb phase (a factor 
« 20), but they do not become worse as /3 increases, as we test in a 12*^ simulation at 
/? = 3.0. 

On the other hand, the geometric algorithm seems to perform much better for the 
imaginary plaquette-plaquette correlation in the confining phase, whereas in the Coulomb 
phase all the advantage vanishes. Again, our 12^ computation at /? = 3.0 reveals that the 
ratio slowly decreases as j3 increases (being f» 0.8 at ft = 3.0). 

Of course, this analysis assumes that both algorithms have no ergodicity problems. We 
must be careful to start from a cold configuration when running the heat-bath simulations 
in the Coulomb phase, in order to avoid metastable states which could spoil the reliability 
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Figure 7: Correlation function of the real part of the plaquette versus plaquette-plaquette distance 
in lattice units, for the four-dimensional U(l) lattice gauge model in the Coulomb phase {(3 = 1.1). 
Beyond distance 4, the error became far larger than the expectation value of the correlation. The 
lattice volume was 16^. 

of the simulation. 

5. The three dimensional Ising gauge model 

Let us finally come to the point of critical slowing down: This is a major issue, as any 
improvement in this field can be of paramount importance in term of the cost of large scale 
simulations of (statistical) systems at a critical point. Beating critical slowing down is one 
of the main motivations in the development of new Monte Carlo algorithms. 

Typically what is found in Monte Carlo simulations of system both in statistical physics 
and gauge theories is that the autocorrelation time r diverges as we approach a critical 
point, usually as a power of the spatial correlation length: r ~ where ^ is the correlation 
length and z is a dynamical critical exponent. Usual local algorithms have values of z 
around 2, maJiing it very inefficient to simulate close to the critical point. For spin systems 
there are well known cluster algorithms with much smaller z. Previously published results 
[19] on an algorithm similar to ours, but applied to a non-gauge model, have claimed a 
similarly smaller vahic for z. Having also this motivation in mind, we have investigated 
the autocorrelation properties of our numerical scheme on the critical point of a system 
that undergoes a second order phase transition (with diverging correlation length). Our 
model of choice has been the three dimensional Ising-gauge model. We have performed 
extensive simulations in the critical zone of this model for several values of the lattice size 
(and hence correlation length) , using both the geometric algorithm and the standard Monte 
Carlo approach, the latter known to have a lower bound for the autocorrelation exponent 
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Figure 8: Correlation function of the imaginary part of the plaquette versus plaquette-plaquette 
distance in lattice units, for the four-dimensional U(l) lattice gauge model in the confining phase 
{/3 = 0.9). Beyond distance 4, the error became far larger than the expectation value of the 
correlation. The lattice volume was 16^. 



z equal to 2, a value typical of all local algorithms. For lattices up to L = 24 we have in 
all cases more than 5 x 10^ Monte-Carlo iterations, which increase to more than 1 X 10^ 
for L = 32, 48, and to more than 4 x 10^ iterations for the largest lattice L = 64. 
For an observable O we define the autocorrelation function p(t) as 

Pit) = iiOii)-0^)iOii + t)-OB)) 

where Oa = {0{i)), Ob = {0{i + 1)), a\ = ((0(i) - Oa?), a| = ((0(i + 1) - Ob? 
and () denotes average over i. We then define the integrated autocorrelation time by 

^ N-f 

r = p(0) + 2^p(i)-r^ (5.2) 

t=i 

where N is fixed, but with N < ?,t and N < 10% of the total sample. In Fig. 11 we 
report the results for the integrated autocorrelation time of the plaquette versus lattice 
size in logarithmic scale for both algorithms. The results of our simulations hint to a 
different asymptotic behaviour of the autocorrelation time, although with our present data 
we cannot obtain a conclusive result. The results for the heat-bath algorithm seem to fall 
nicely on a straight line, which would correspond to a simple exponential dependence of 
r on L, with z = 2.67 it 0.08, but the geometric algorithm presents a more complicated 
behaviour, as well as larger errors. There are signs that the asymptotic behaviour might 
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Figure 9: Correlation function of the imaginary part of the plaquette versus plaquette-plaquette 
distance in lattice units, for the four-dimensional U(l) lattice gauge model in the Coulomb phase 
(/3 = 1.1). Notice the different behaviour of the algorithms at large distances. The lattice volume 
was 16*. 

be better than for the heat-bath, but much more extensive simulations, outside the scope 
of this work, would be needed to get a definite value for z. 

6. Conclusions and Outlook 

Our main motivation is the sign problem, particularly in QCD at finite chemical potential, 
biit which also appears in other systems of interest, for example when a Q vacuum term is 
present. New ideas are clearly needed in order to make significant advances in this problem, 
and one possibility is the development of new simulation algorithms that might circumvent 
the difficulties of conventional approaches. 

We have developed a geometric algorithm, based on the strong coupling expansion 
of the partition function, which can be applied to abclian pure gauge models. We have 
checked in the U (1) model in 3 and 4 dimensions that the algorithm can be implemented 
efficiently, and is comparable with a standard heat-bath algorithm for those models. It 
seems however that the geometric algorithm does not suffer lack of ergodicity due to the 
presence of vortices, as can be the case for heat-bath, depending on the starting point. 

We have also studied the algorithm in the 3 dimensional Ising gauge model at the 
critical point, where we have seen hints that the asymptotic behaviour of the geometric 
algorithm may be better than standard heat-bath. This would be very interesting, because 
in contrast to spin systems, where there exists cluster algorithms that can greatly reduce 
critical slowing-down, to our knowledge no similar algorithm is known for gauge systems. 
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Figure 10: Ratios of figures of merit for different observables between geometric and heat-bath 
algorithms. Re Corr/{=2 and Im Corrfl=2 stand for Real and Imaginary plaquette-plaquette corre- 
lations at distance 2. 
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Figure 11: Autocorrelation times at the critical point (of each algorithm) versus lattice length; 

boxes stand for standard algorithm results, with a linear fit to guide the eye, while circles represent 
the results of the geometric algorithm. The errors were obtained by a jack-knife procedure. 



Our results are however not enough to establish this, and much more extensive simulations 
should be done to clarify this point. 



The algorithm can be extended to include fermions, and this constitutes work in 
progress. In this case there is a sign problem, and one question we want to answer is 
whether such problem is severe or mild. 
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